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Abstract. Quantum mechanical entanglement can exist in noisy open quantum 
systems at high temperature. A simple mechanism, where system particles 
are randomly reset to some standard initial state, can counteract the 
deteriorating effect of decoherence, resulting in an entangled steady state far 
from thermodynamical equilibrium. We present models for both gas-type systems 
and for strongly coupled systems. We point out in which way the entanglement 
resulting from such a reset mechanism is different from the entanglement that one 
can find in thermal states. We develop master equations to describe the system 
and its interaction with an environment, study toy models with two particles 
(qubits), where the master equation can often be solved analytically, and finally 
examine larger systems with possibly fluctuating particle numbers. We find that 
in gas-type systems, the reset mechanism can produce an entangled steady state 
for an arbitrary temperature of the environment, while this is not true in strongly 
coupled systems. But even then, the temperature range where one can find 
entangled steady states is typically much higher with the reset mechanism. 



PACS numbers: 3.65.Yz, 03.67.Mn 

1. Introduction 

In quantum information theory entanglement between parts of a system has 
been identified as the key resource that can possibly make quantum information 
processing more powerful than classical information processing. Entanglement can 
also be a resource for long-distance quantum communication or distributed quantum 
computation, and it is at the heart of some quantum communication protocols. But 
entanglement is fragile under the influence of environment induced decoherence. All 
engineering hence thrives to better control and manipulate the quantum information 
stored in the system while keeping the detrimental effects of decoherence low. 

In nature, on the other hand, we mostly find less controllable systems, especially 
if the system size becomes macroscopic as in gases, fluids, solids or even biological 
systems. Since these systems are usually open, noisy systems at possibly high 
temperatures one expects that environment induced decoherence will erase all 
entanglement between system degrees of freedom. This reasoning is true except for 
three cases. 

First, the environment and its coupling to the system could be special in a way 
that it creates rather than destroys entanglement. However, it is unlikely to find such 
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an environment in nature where usually thermalization dominates, and we will only 
briefly touch upon the subject of such environments in this paper. 

Second, if the system has an entangled ground state, as many solid state systems 
do, its thermal state will be entangled in a certain temperature range above zero by a 
continuity argument. Coupling to a heat bath drives a system into its thermal state. 
But there is a temperature threshold for the bath above which the thermal state will 
be unentangled. 

Third, the system might have a built-in entropy drain, meaning that the 
correlations with the environment are, by one way or another, erased such that the 
system can re-build entanglement through its quantum mechanical interactions. This 
entropy drain may even be local to exclude the trivial cases where entanglement is 
simply "pumped" into the system, e.g., by injecting fresh, entangled Bell pairs. 

In [I] we proposed such a local entropy drain in form of a reset mechanism, where 
system particles are randomly replaced by particles in some standard, mixed state 
of sufficiently low entropy. Note that such a mechanism cannot create entanglement, 
on the contrary, it erases any entanglement that might still be present between the 
particle that is reset and the rest of the system. Only the interplay with the system 
Hamiltonian can lead to entanglement in a steady state that is possibly far from 
thermodynamic equilibrium. This reset mechanism was studied for a toy model with 
two qubits, where analytic solutions could be obtained. Also a multipartite scenario 
for a (simplified) gas model was discussed, and further generalizations were suggested. 
By gas-type systems we mean systems in which the decoherence processes act locally 
on the system particles, by strongly coupled systems we mean those where the 
decoherence processes act globally. To be more precise, local decoherence processes are 
those, which induce transitions between the eigenstates of the local, free Hamiltonian 
alone, while global decoherence processes induce transitions between eigenstates of the 
total Hamiltonian. 

In this paper we review the key idea of a reset mechanism but provide more 
in-depth material than in [1]. We elaborate on the generalizations suggested in [T], 
namely on the influence of local entropy drains on the dynamics of entanglement and 
on the steady-state entanglement in gas-type systems as well as in strongly coupled 
systems. 

We prove that the master equation describing the evolution of the system coupled 
to a heat bath and subject to a reset mechanism is of Lindblad form and hence 
generates a completely positive, i.e., physical map. We analytically solve the master 
equation for small systems of two spins with special interaction Hamiltonians, which 
enables us to illustrate the main features of the reset mechanism. In particular we 
show the following. 

(i) Steady-state entanglement in systems with reset mechanism is different from the 
entanglement in thermal states. 

(ii) In strongly coupled systems with constant coupling steady-state entanglement 
with reset can exist for higher temperatures than the entanglement in the thermal 
state, which is the steady state without reset. 

(iii) In gas-type systems steady-state entanglement with reset can exist even for 
arbitrary temperatures. 

These features are not due to the specially chosen interaction Hamiltonians 
and decoherence processes. Wc demonstrate that the above properties are almost 
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independent of both. One can also relax the conditions on the reset states and 
take mixed states with sufficiently low entropy instead of pure states. Finally, a 
generalization to larger system sizes, possibly even with fluctuating particle numbers, 
still leads to similar results. Hence, the reset mechanism is at the same time simple 
and generic. 

We remark that in cavity QED an incoherent generation of entanglement has been 
proposed, which bears resemblance to our work [2]. There, an atom couples to two 
leaky optical cavities and is driven by a white noise field. This incoherent driving can, 
when the atom is finally traced out, result in entanglement between the cavity modes. 
Entanglement is generated for intermediate cavity damping rates and intensities of 
the noise field, an effect labeled "stochastic resonance" in [2 [3]. We believe that 
this effect is more correctly interpreted as an example for a reset mechanism. In a 
subsequent work [4] strongly related to [2], one single cavity entangles two atoms, 
giving yet another example for a reset mechanism even closer to the setups of this 
paper. 

The reset mechanism is certainly not a preferred way to actively protect 
entanglement and mostly cannot even be compared to such strategies, but, because of 
its simplicity and generality, there is hope that such a mechanism may ultimately 
be identified in natural processes leading to an increased understanding whether 
entanglement can play a role in systems at high temperatures. 

The paper is organized as follows. We first concentrate on simple models with only 
two particles (qubits). In section [2] we motivate the description by a master equation, 
explain in which cases the model is valid, and study several specific Hamiltonians 
and noise channels analytically and others numerically. We also compare entangled 
steady states resulting from a reset mechanism to entangled steady states resulting 
from special choices of interaction Hamiltonian and decoherence process. We show in 
section [3] that we can find the same features in strongly coupled systems, and we give 
the conditions to be met by the reset mechanism such that entangled steady states can 
exist. Then, in section [4j we extend the model to include more qubits and discuss the 
meaning of different kinds of entanglement that we use. Finally we give a summary 
of the results in section 

2. Gas-type systems 

In this chapter we discuss a toy model with only two particles, which we take as spin- 
1/2 systems or qubits for simplicity. The toy model shows all the features that we will 
later find in larger systems and it has the advantage that we can show many results 
analytically leading to an increased understanding of the involved processes. We will 
formulate the equations for an arbitrary number N of qubits, so that we can refer to 
them later in section [4] 

2.1. Master equation for gas-type systems 

In a gas particles are weakly coupled in the sense that most of the time they do 
not considerably interact with each other unless they collide. In the meantime they 
only feel their local, free Hamiltonian and are subject to individual, local decoherence 
processes, e.g. through interactions with thermal photons (radiative damping). If 
we pick a subset consisting of N gas particles and consider these as the system, 
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collisions with the remaining gas particles are another source of decoherence (non- 
radiative damping or dephasing). In a master equation that models this gas- type 
system we replace the original, time-dependent collision Hamiltonian by an averaged, 
time-independent interaction Hamiltonian. Since the interaction Hamiltonian does not 
modify the energy landscape in this model, the local, radiative decoherence processes 
tend to drive the system to the thermal state of the free Hamiltonian, for which we 
choose the form 

N 

tf frcc ^/2]>>i 4 ). (1) 

i=l 

We leave the interaction Hamiltonian H unspecified for the moment. For two qubits, 
we will often use the Ising Hamiltonian 

tfi S i„g = g4 1] 4 2) (2) 

for analytic discussions, whereas more complicated Hamiltonians will be treated 
numerically. We write the total Hamiltonian as H to tai = H + H[ ICC such that the 
master equation is 

P = ~i [-Htotai, P] + Aioisc/O, (3) 

where £ n oisc is a Liouville operator representing the noise channels. We describe the 
noise channels by the Lindblad operator [5J [5] 

N 



D 



^ 2 "(1 - s)\a\'al_'p + pa\>a'_> - 2a y _'pa y + '\ 
i=l 

-—s[o- y _'o-yp + pal J o-y - 2o-ypal J ] [p-a\ } pa\ } \ (4) 

where a± — (cr x ± ia y )/2 and the as are Pauli operators. Parameters B and C give 
the decay rate of inversion ( n" ) and polarization (<J±) under the action of C no [ se , 
and s = lim t ^ oc ((l + <r z )/2) t — (e w/3 + G [0, 1] depends on temperature, where 
s = 1/2 corresponds toT = l//3 = oo (we set the Boltzmann and Planck constant 
equal to one). The definition of s stems from laser physics where inversion occurs 
corresponding to "negative temperatures" . Many authors use n = l/(e w/3 -l) instead 
of s and n+1 instead of (1— s). Then, no negative temperatures are possible. The noise 
channel is derived assuming certain approximations, e.g. the Markov approximation. 
Note, however, that this is not an essential assumption as we will demonstrate later 
in an example (see Figure [2] and related text). 

An important special case of H]), obtained by setting B — and C — 2j, is the 
dephasing channel 

JV 

£dephP = 7E[ CT *V^-p], (5) 
i=l 

well-known especially in its integrated form as a completely positive map £ ^i)deph(p) = 
pp+ (l—p)/2 (a^ po~^ + p^j with p = exp{— 2"ft}. As with the Ising Hamiltonian @, 
we will often use the dephasing channel for analytic discussions because of its 
simplicity. 

We are interested in the steady state of this master equation for N = 2 qubits 
and the question whether there is entanglement in this state. At this point we simply 
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state the following results since we will later solve a more general master equation 
that contains equation as a special case. 

As an easy example we start with the Ising interaction Hamiltonian The 
steady state will be the tensor product of the thermal states of each free Hamiltonian 
iTfree = cj/2 <j£ since this state commutes with i?i s i ng and since dephasing noise does 
not change the diagonal elements of the density matrix. In conclusion, we have the 
unentangled steady state 

Pstoady = diag(s 2 , s(l - s), s(l - s), (1 - s) 2 ). 

When can we hope to find an entangled steady state? As we know the radiative 
decoherence processes drive the system into the thermal states of the free Hamiltonian. 
If the interaction Hamiltonian can entangle these states we may find an entangled 
steady state at least for low temperatures. As an example, consider an interaction 
Hamiltonian H — gcri^a^ while the other terms stay the same as in the example 
above. We set C — 1/2 B (no non-dissipative processes). The steady-state density 
matrix is then 

/ g 2 +s 2 (B 2 +4uj 2 ) g(2s-l)( l B+2uj) 
B 2 +4{g 2 +u 2 ) U U B 2 +4(g 2 +ui 2 ) 

U B 2 +4(g 2 +u 2 ) U ' 







g 2 -(s-l)s(B 2 +4u J 2 ) 
B 2 +4(g 2 +u> 2 ) 

g(2s-l)(2w-iB) n n g 2 +B 2 (s-l) 2 +4(s-l) 2 ^ 

\ B 2 +4(g 2 +J 2 ) U U B 2 +4(g 2 +uj 2 ) / 

This density matrix can be entangled. We measure the entanglement between 
two sets of qubits by the negativity pj [5] , which is given with respect to a bipartition 
A- A as Ma = (IIp Ta I|i — l)/2, where Ta means the partial transpose with respect 
to A. For two qubits, we omit the label A since there is only one bipartition, and 
the negativity can assume values between zero (separable state) and 1/2 (maximally 
entangled state). The reason why we choose the negativity as a measure throughout 
the paper is that we will use a generalization thereof in the multipartite case where 
the generalization of other entanglement measures might be hard to compute. 

The negativity of the state above is 

N = Max{0, [ (B 2 + 4 (g 2 + lo 2 )) ({s - l)s (B 2 + 4w 2 ) - g 2 ) 

+g\l 2s\ (B 2 + 4c 2 ) 1/2 (B 2 + 4 (g 2 + c 2 )) ] (B 2 + 4 (g 2 + u 2 )y 2 }, (6) 

which can be larger than zero but will always vanish for high temperatures of the bath, 
s — > 1/2 (see Figure [1]). We see that the steady state of the master equation © can 
be entangled for specially chosen interaction Hamiltonians, but only below a certain 
temperature threshold. Furthermore, if the free, local Hamiltonian is too strong, i.e., u> 
dominates by far all other parameters, there is also no entanglement. This statement 
applies to other models involving the Hamiltonian |T]) as well, and, accordingly, uj 
should have the same order of magnitude as the other parameters. 



2. 2. Example for a gas-type system: the spin gas 

Spin gases [9l [10] are an example for such gas-type systems. A spin gas is a system 
of quantum spins with stochastic, time-dependent interactions. A physical model 
of a spin gas is a system of classically moving particles with additional, internal 






Figure 1. Entanglement as measured by the negativity in the steady state of 
equation J3|l with interaction Hamiltonian H = gpjb^tr^ and free Hamiltonian 

Hfree = w/2 JDiLi o* • The decoherence processes are given by equation Kl . 
We choose the decay rate of inversion, _B, as inverse unit timescale and the other 
parameters as C = B/2 (no dephasing noise) and u = B. The parameter g on 
one of the axes is also measured in units of B, whereas s is dimensionless. 

spin degrees of freedom. Upon collision, these quantum degrees of freedom interact 
according to some specified Hamiltonian. In [9, 10] the interaction Hamiltonians were 
chosen locally unitarily equivalent to the Ising interaction leading to a description 
in terms of weighted graph states. Hence, in such spin gases, classical kinematics 
drives the evolution of the quantum state, and also the decoherence of arbitrary 
probe systems put into the gas and subjected to interactions with it. In general, 
multiple non-consecutive collisions of particles are possible. The spin gas remembers 
its whole interaction history, and it provides a microscopic model with non-Markovian 
decoherence. 

Assume that we have two selected gas particles (e.g. another species) that we con- 
sider as the system, while the other gas particles act as the environment that induces 
decoherence when interacting with the system particles. The rare interactions between 
the two system particles happen only during the short times when they collide. In the 
longer times in between they are not coupled and subject to local decoherence pro- 
cesses, i.e., interactions with the environment gas particles. The induced decoherence 
processes are equivalent to dephasing channels (corresponding to B = 0, C = 27 in 
equation In such a situation, any entanglement between the two qubits that may 
either have been present initially or have built up on a short time scale will eventually 
be destroyed by the interactions with the other gas particles [TJ. 



The spin gas with reset mechanism 



Entanglement and its dynamics in open, dissipative systems 



7 



For the moment, we stick to the toy-model with only to selected gas particles. 
Imagine now that the two particles can, at a certain rate, leave the box in which the 
gas is contained and are instantly replaced by fresh qubits that are in a standard mixed 
state with sufficiently low entropy. Instead of a replacement of system particles one 
can equivalently picture a measurement of the particle and a subsequent preparation 
in this standard state. Note that the last step need not be an active procedure but 
can, e.g., result from a spontaneous decay to this state. We call both procedures a 
reset mechanism. Certainly, by a reset, we did not introduce entanglement into the 
system always consisting of two qubits. On the contrary, any entanglement that might 
have been present between the particle that has left the box and the one that is still 
inside leads to a description of the latter by a more mixed density matrix (closer to the 
identity). But the advantage is that we have lowered the local entropy of the system 
since the new particle has no correlations with the environment. This new particle can 
then become entangled with the other one on short time scales. We said above that in 
a spin gas with zero rate of qubit exchange the steady state will not be entangled. For 
infinite exchange rate the system would always be in a pure (or a standard separable) 
state and there is also no entanglement. If, however, the rate at which the qubits leave 
the spin gas is in a certain intermediate parameter regime one can hope that there 
is entanglement in the system on average. Here, averaging means taking the mean 
density matrix of many simulation runs. Later, the solutions to master equations are 
assumed to resemble the evolution of such a mean density matrix and, for explanations 
of certain (entanglement) features, this picture will sometimes be invoked. Note that 
there are also other ansatzes. In [11] the solution of the master equation represents a 
smoothed version of a single simulation run, where smoothing is achieved by a time- 
integration kernel. The solution of the master equation does then not follow the rapid 
changes of the single density matrix, but sees only the slower changes resolved by a 
so-called coarse-grained timescale, which is related to the support of the integration 
kernel. Figure [2] indeed shows entanglement in steady states in a simulation of a spin 
lattice gas with an Ising-type interaction [9] [10] . 

The above scenario with two qubits might seem a little artificial. However, if we 
extend the setup to systems with more qubits and allow fluctuating particle numbers 
we can drop the requirement that selected particles must be instantaneously replaced. 
We can regard the spin-gas particles that do not belong to the system as a "hot bath" 
introducing decoherence and destroying entanglement in the system, and the source 
of qubits in standard states as a "cold bath" that can counteract the effect of the hot 
bath and preserve entanglement in the steady state. The analogy with a cold bath 
has some limitations as we will point out in section [3] We will deal with multipartite 
qubit systems in section [4] For the moment, we will stick to the two qubit system, for 
which we can find a master equation that reflects the properties of the "cold bath" and 
that we can solve analytically in certain cases. Observe that the master equation will 
again incorporate the Markov-assumption, whereas the spin gases are non-Markovian 
systems [9] [10] and also partly have non-local decoherence processes. The essential 
features on the other hand will be qualitatively the same in some parameter regimes 
of the spin gas, where these effects play a minor role. 

In the next subsection, our goal is to transfer the idea of a reset mechanism from 
the specific example of a spin gas to a description in form of a master equation, suitable 
for any spin system. 
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Figure 2. Steady-state entanglement between two selected qubits as measured by 
the negativity in a spin lattice gas (8x8 lattice, 2+18 qubits). The probability that 
a particle is exchanged for a fresh one in one time step of the simulation is plotted 
on the horizontal axes. Hence, the value 1 corresponds to an infinite exchange 
rate. The special qubits interact 1000 times stronger with each other than with 
the 18 qubits that form the environment, i.e., as physical particles, they are e.g. 
of a different kind than the environment particles. The density matrices from 
which the negativity is derived were averaged over 10000 simulation runs. Other 
parameters: initial distance between special qubits: 1, interaction phase picked 
up during a collision between them: iji = 0.1; interaction phase for interactions 
with environment spins: <f> = 0.0001; for details see |91|10|. 



2.3. The master equation with reset mechanism 

Compared to ([3]), the master equation that models a gas- type system with reset 
mechanism has an additional term £ reset , which we describe as follows. With some 
probability rSt particle i, i = 1...N, is reset during the time interval St to some specific 
state |xi)- The other qubits are left in the state trip. The change in the density matrix 
during the time St due to C reset is (Sp) reset = rStY, l= i^(\Xi){Xi\^iP- P) = StC reset p. 
Observe that the time interval must be longer than the timescale of any of the involved 
processes but short enough so that we can replace it by the time differential to obtain, 
for the rate of change p — dp/dt, the following master equation: 

N 

P = -i[H t otal,p] + CnoiseP + ^2r(\Xi)(Xi\tiiP " P) (?) 

i=l 

Before we proceed to discuss the solution of ([7]) let us establish that the problem 
is well-defined, i.e., that the master equation leads to a completely positive map, 
which is true when the master equation is of Lindblad form. For the noise part 
this is known, so we have to bother only about the reset part. Since C rese t = 
X}j=i r {\Xj)(Xj\^ r jP — p) is focal we have to show that each summand is of the form 

Sm,n=i 4n(kraV> °n \ + [ cr m ) , P^P] ) where the er (j) s are Pauli operators and U 
must be positive (semidefmite) matrices. We expand p and |Xj)(Xj| m t ne cr-basis 

as P = J2l u ...,k N =o Ofe,..,**^ ■ ■ ■ a Z } and \Xj)(Xj\ = Eq=o H^- We insert these 
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expressions into C reset and also into the Lindblad-expression, collect the coefficients 
that belong to each cr-matrix of p using the scalar product, and compare the coefficients 
dfci in each expression, which leads us to a simple linear system of equations for the 
L J mn . Solving this system of linear equations we obtain 

(I -H l M \ 

L J =r[ \iV 3 I -\ib{ (8) 

The eigenvalues of U arer/8,r/8 (l±2((b\) 2 + (b J 2 ) 2 + (fr^) 2 ) 1 ^. Since we assumed 

\Xj)(Xj\ to be pure, we have (6{) 2 + (b 3 2 ) 2 + (b^) 2 — 1/4 and eigenvalues r/8,0, r/4. 
We note that also a mixed reset state would be fine to ensure that the L 3 are positive 
semidefinite. Because the sum of positive matrices is a positive matrix, and because 
we know that the noise terms also have positive L-matrices, we have shown that the 
master equation is of Lindblad form and preserves the positivity of the density matrix. 

Up to this point, we have modeled interacting, gas- type systems coupled to a 
noisy environment. We have described a toy model consisting of only two particles 
by a master equation and compared predictions about the entanglement properties of 
steady states to simulations with a spin gas as an example for such gas-type systems. 
We have seen that in general there will be no entanglement in the steady state. We 
have extended the example of the spin gas by allowing particle exchange with a "cold 
bath" of particles in standard states (or an equivalent reset mechanism), and we have 
found that steady states of such systems can be entangled. We have derived a master 
equation that models systems with reset mechanism and have proved that the master 
equation is of Lindblad form. In the following subsection we study the solutions of (JT]) . 

2.4- Solution of the master equation for the gas-type model 

In principle, the solution to the master equation ([7]) with noise channels as in 
equation (j4|) is simple. The equation is of the form 

p = Cp 

with solution 

p(t)=e ct p(0). 

Mapping p to a column vector C containing the 16 coefficients Cbooo, Coooi, • ■ ■ , Cuu 
of the density matrix and accordingly mapping the Liouville operator C to a 16 X 16- 
matrix A we get the equivalent 16 coupled linear differential equations C — AC with 
solution C{t) — e At C(0). To compute the matrix exponential we need the spectral 
representation of A, i.e., we must solve the eigenproblem AC\ = XC\. Observe that 
the steady state (if it exists) is given by the eigenvector C corresponding to the 
eigenvalue A = 0. 

In the following we will first analyze the solutions for an Ising interaction 
Hamiltonian and later generalize to generic cases. 

2.4-1- Ising Hamiltonian. We specialize to the Ising Hamiltonian ([2]) as (effective) 
interaction Hamiltonian, H = Hi sms , and to a specific reset state, namely |+)(+| 
for both qubits. The free Hamiltonian is £/f rec = w/2 YnLi as before. We 
can solve the problem through spectral decomposition of the Liouville operator 
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C(p) = —i[H tota \,p] + AioiseP + A-csct/ ! but the expression for the corresponding 
matrix e At is very lengthy. 

One obtains shorter expressions if one does not solve all 16 differential equations 
at once through the matrix exponential, but step by step, since not all differential 
equations are coupled. Still, we have chosen to move the solution derived in this way 
to | Appendix A| not to overburden the text with technical details. 

For illustration, we will restrict the noise to the special case of a dephasing chan- 
nel (0 in the following. 

Dephasing channel 

As pointed out above, the solution is given by the spectrum of the total Liouville 
super-operator defined by p — Lp and its corresponding eigenvectors. In the case of a 
dephasing channel one obtains the eigenvalues 

{0, -r, -2r, -2(r + 2 7 ), -2(r + 2 7 ± iu), 

-l/2(3r + 4 7 + y / -l6g 2 + r 2 ± 2ioj, -l/2(3r + 4 7 - ^-16g 2 + r 2 ± 2iu} 

with multiplicities {1, 2, 1, 2, 1 + 1, 2 + 2, 2 + 2}, respectively. The eigenvector belonging 
to the eigenvalue represents the density matrix in the steady state, and we will 
come back to this matrix in the next subsection. The full solution, derived by 
solving the differential equations in a step-by-step manner as explained above, is given 
in |Appendix A| 

To demonstrate the time-evolution of an initial density matrix governed by the 
master equation we plot the entanglement between the two qubits as a function of 
time. Note that Figure [3] is based on the analytic solution given in |Appendix A[ 
and the plot shows how the negativities of different initial states approach the final 
negativity of the steady state. We choose 7 -1 as unit timescale (setting 7 = 1). 
The parameters r, <?, lo here have the fixed values IO7, 5 7 , 57, respectively. The 
initial states are weighted graph states [91H2] with density matrix U((p)\+) (+\W(<p) 
where U(ip) = diag(l, 1, 1, e lip ). Through the parameter <p we can continuously tune 
the entanglement in the initial state from the product state |+)(+| for ip = to the 
maximally entangled, Bell-equivalent state for ip — it. States that are initially highly 
entangled are first driven into separable states before the steady-state entanglement 
value is approached from below. Vice versa, an initial product state gets highly 
entangled first, before the steady-state value for the entanglement is reached from 
above. As we said earlier, to display the full analytic solution of ([7|) for more general 
noise channels would be quite space-consuming. We will not present it since we 
are primarily interested in the entanglement properties of the steady state. In the 
following we will discuss these properties for the master equation with general local 
noise channels. 

2.1^.2. General steady states with Ising Hamiltonian. As we have seen, any initial 
state of the density matrix evolves exponentially fast into a steady state on a 
characteristic timescale given by the largest non-zero characteristic exponent (or the 
smallest in absolute values since they are negative) . The characteristic timescale thus 
depends on the parameters of the master equation, too. The steady state is also a 
function of these parameters. 

To smooth the presentation, we have again transferred the steady-state solution 
of j7|) with Ising Hamiltonian, local noise channels as in ([4]) and reset states = |+) 
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Figure 3. (a) Time development of entanglement for different initial states. The 
2-axis displays the entanglement as measured by the negativity. The initial states 
are weighted graph states, characterized by a parameter ip (see text). The unit 
timescale is 7 — 1 , t is measured on this timescale, and the parameters of the master 
equation are r = IO7, g = 5-f, u> = 57. (b) Cut through the 3D plot at tp = 
(orange curve) and (p = ir (blue curve). 



to I Appendix B| Here, we illustrate the solution with a plot. 

In Figure H] we choose B~ x as unit timescale and fix the values C — B, ui = 20B 
and the dimensionless parameter s = 0.1. Then, for certain values of r and g, 
measured in units of B, we see entanglement as measured by the negativity in the 
steady state. To get entangled steady states when g becomes small, we have to go 
to higher reset rates r. However, g cannot be arbitrarily small. There is a weak 
coupling threshold below which no reset rate can ensure entanglement in the steady 
state. This threshold depends on the parameters of the decoherence processes, and 
its existence is intuitively clear. If the decoherence processes simply dominate the 
entangling processes, then no entanglement can be created by any means. To see 
this better, we momentarily put uj = for simplification. Let us also turn our 
attention once more to the dephasing channel as a special case of the noise terms 
of the quantum-optical master equation. Then, in the steady state, the anti-diagonal 
coefficients are all the same, namely 4( r+2 7)(2g 2 +(r+7)(r+27)) ' ^ nc other off-diagonal 

elements are C 00 oi = C0010 = Cg m = C* 011 = ^2^f^^+2jT) ' and the dia S onal 
elements all have the values j . All other matrix elements are given by the Hermiticity 
of the density matrix. We compute from the above expressions for the density matrix 
the following analytic expression for the negativity in terms of the parameters g 
(Hamiltonian interaction), 7 (strength of the dephasing channel), and r (reset rate): 

2 7 ( r + 7 ) 2 + g 2 (r + 2 7 ) - r(r + 2-/)g 



M = max{0, 



-} 



(9) 



2(r + 2 7 )[2 ff 2 + (r + 7 )(r + 2 7 )] 

Equation ([9]) contains the full information about the entanglement properties of the 
two qubits. Note that J\f — AfQj, r) depends, in fact, only on two parameters, g = g/ 7 
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Figure 4. Entanglement in the steady state of the master equation JZ}. The 
entanglement is measured by the negativity, the unit timescale is B —1 , the 
parameters r and g (in units of B) are on the axes, while the other parameters 
are C = B, lu = 20B, and s = 0.1. 

and f = r/7. 

In Fig. Ela) we see a plot of the negativity function N. The key feature is the 
color-coded region in the r-g-plane with steady-state entanglement, where a darker 
color indicates higher entanglement. The entangled region is bounded by the red line 
given by one of the roots of the non-trivial part of equation ([9|). Outside of this region, 
the state is separable (white area). The entangled region approaches asymptotically 
the straights g = 2-f and g = r plotted black in Fig. 02a). The asymptotic line 
g = 27 is independent of r and simply tells us that, in the weak coupling regime, 
decoherence/noise will always triumph over the Hamiltonian part that tries to create 
entanglement as pointed out before. That is, as a necessary condition, we need to be 
above this threshold to observe entanglement. Three lines are marked in the colored, 
entangled region: 

(i) The upper, white line is the maximum in (/-direction (at constant r). 

(ii) The lower, white line is the maximum in r-direction (at constant g). 

(iii) The middle line in black is the straight g = r/(l + \/3). To this middle line the 
upper and lower white curves go asymptotically for large g, r. 

The global maximum of the negativity is on this middle line at infinity with a value of 
approximately 0.0915, about 20% of the maximally possible value. The darkest, most 
entangled area in our plot has negativity approximately 0.068. Fig. EJb) shows a cut 
at g = 57 through the color-plot. Most notable is the existence of a threshold value 
for r/7 above which entanglement is present in the steady state. 

2.4-3. Steady state entanglement with and without reset. In this part we want to 
compare steady-state entanglement that is due to the reset to entanglement that is 
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Figure 5. (a) Separable states (white area) and entangled states (colored area) 
in the r-g-plane, where r is the rate of the reset process, and g is the coupling 
strength in the Ising Hamiltonian, and we use 7 _1 = 1 as unit timescale. The 
color encodes the amount of entanglement measured by the negativity: the darker 
the area, the more entanglement is present. For a discussion of the other lines, 
please see the text, (b) The second plot shows a cut for constant g = 5-f. 



due to a special combination of interaction Hamiltonian and decoherence process as 
in subsection [2TTJ There, the interaction Hamiltonian is 

while the free Hamiltonian are given by (TT]) . The decoherence process (Q| is determined 
by the parameter choice C = B/2. In l2.11 it has been established that the steady state 
is entangled, but only in a finite temperature range. 

If we add the reset mechanism to the master equation of this example we can 
show that the steady state is entangled for arbitrary temperatures. As reset states we 
choose the eigenstate |1) = — a z \l) of the Pauli operator <r z . 

The steady-state density matrix of the master equation has now the matrix 
elements 

B 2 s 2 uj 2 + (B + 2r) {(B + r)g 2 + B 2 (B + 2r)s 2 ) 



Ci 



0000 



(B + r) {{B + r)u 2 + (B + 2r) (4g 2 + (B + r)(B + 2r))) 
Coioi = Cioio = {(B + r) ({B + r)u 2 + (B + 2r) (4g 2 + (B + r)(B + 2r)))} _1 
x {(B + 2r) ((B + r)g 2 - B 2 (B + 2r)s 2 + B(B + r){B + 2r)s) - Bs(sB - B - r)uj 2 } 

Cim = {( B + r ) (i B + r V 2 + (B + 2r) (4 5 2 + [B + r){B + 2r))) } _1 
{(-sB + B + r) 2 uj 2 

+(B + 2r) [B 2 (B + 2r)s 2 - 2B{B + r)(B + 2r)s + (B + r) (g 2 + (B + r)(B + 2r))] } 
g(2sB-B-r)(i(B + 2r)+u>) 



0011 



(B + r)cu 2 + (B + 2r) {Ag 2 + {B + r){B + 2r)) 

while all other coefficients are zero or are given by Hermiticity. The negativity of this 
density matrix is 



1 



((B + 2r) 2 +uj 2 ) (-2sB + B 



Af max^O, 4 I - + ^ ^ + ^ + + ^ + + + ^ 
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g 2 (-2sB + B + r) 2 ((B + 2r) 2 



((B + r)uo 2 + {B + 2r) (4g 2 + (B + r)(B + 2r))y 



1 



(11) 



We plot the function (fTTjl in Figure [6l see figure caption for details. We observe 
two features: 

(i) At r = we are back to the situation of subsection 12.11 where entanglement 
vanishes above some temperature threshold (remember: T — > oo as s — > 1/2). 

(ii) There is a threshold reset rate r, above which an entangled steady state exists for 
arbitrary temperature. 

From the coefficients of the steady-state density matrix (fTU|) we also see that the 
entanglement created by the reset stems from a different density matrix than the 
entanglement present without reset. In Figure [5] this is visible in the region of small 
r, where the reset tends to destroy this latter entanglement. Then, for larger r, the 
effect of the reset mechanism kicks in. 




Figure 6. (a) Negativity for 2-qubit system with H tota i = g(Tx u<K x^ +w/2(<rj ^ + 
cri 2 ^) as a function of the reset rate r and temperature-dependent parameter 
s. The noise is described by the quantum-optical master equation with 
B = 2C = 1. With B^ 1 as unit timescale the other parameters are given by 
g = IB and lu = 2B. The entanglement due to the reset mechanism exists for all 
temperatures, while the entanglement without reset mechanism (r = 0) vanishes 
above a certain temperature threshold, (b) Cut at constant s = (orange curve, 
corresponding to zero temperature) and s = 0.5 (blue curve, corresponding to 
infinite temperature). 



Up to now we studied rather special interaction Hamiltonians. In the following 
we demonstrate the genericity of entanglement that is present in a steady state due 
to a reset mechanism. 



2.4-4- Generic cases. In the previous example we have pointed out that 
entanglement, if present at all without reset, stems from a special combination of 
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Hamiltonian and decoherence process. One may ask whether adding a reset mechanism 
with special reset states is not just as artificial as the choice of special combinations 
of Hamiltonians and decoherence processes. We are now going to show that one 
and the same reset with fixed reset states can lead to steady-state entanglement 
for many combinations of Hamiltonians and decoherence processes. The reasoning 
was the opposite without reset, where only very few combinations of Hamiltonians 
and decoherence processes lead to steady-state entanglement. We can also relax the 
condition that the reset states are pure states to a certain extent. Hence, we show 
that we have found generic features by generalizing the system in various directions. 

(i) The qualitative behavior of the two-qubit model does not depend on the 
particular choice of the interaction Hamiltonian or details of the decoherence model 
other than its local action on individual qubits. Figure [7] shows e.g. steady-state 
entanglement for an XYZ Hamiltonian as function of reset rate r, and decoherence 
described by the noise operator £ no i se . The qualitative behavior is similar to 
Figure (Ub) or Figure and we observe steady-state entanglement even for infinite 
temperature of the bath. 




Figure 7. (a) Negativity for two-qubit system with XYZ interaction and 
magnetic field, H = g{0.7cr ( x 1) + .3cr ( y V ' ' a ( y ] + a^' ' cr^ + Q.5((t^ + u^)) + 

oj/2(az ' + <Jz ))i as a function of the reset rate r and temperature-related 
parameter s. The noise is described by the quantum-optical master equation 
channel {f4} with B = 2C = 1 (B~ 1 as unit timescale). The Hamiltonian 
parameters are g = 2.5B, w = 4B. (b) Cut through the plot. The upper curve 
corresponds to zero temperature (s = 0), the lower one to infinite temperature 
(s = 1/2) of the bath. Curves for any finite temperature lie in between. 

(ii) The idealized reset mechanism we consider can be replaced by a more realistic 
imperfect reset mechanism. In this case, fresh particles are in mixed states with 
sufficiently low entropy rather than in pure states (with entropy 0). Still, the steady 
state turns out to be entangled. When we vary r there is a new, third threshold. 
First, for very small r, there can be an entangled steady state, which is not due to 
the reset and which is present only for a finite temperature range above zero. Second, 
there is one threshold value for r above which the steady state is entangled due to the 
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reset mechanism (for arbitrary temperature). Third, whereas for pure reset states the 
entanglement goes down to zero again only in the limit r — > oo (permanent projection 
to the product of the reset states) , the entanglement goes down to zero for some finite 
r if the reset states are mixed states. This behavior is easy to understand, since it is 
more difficult for the interaction Hamiltonian to create entangled states from mixed 
reset states. The higher the entropies of the reset states are, the smaller is the range 
of the reset rate r for which there is entanglement in the steady state. This range 
can also become zero, so we must demand reset states of sufficiently low entropy. 
Figure ijHJ) clearly shows this new threshold appearing for large reset rates. 




Figure 8. Negativity for two-qubit system with X Y Z interaction and parameters 
as in Figure (0. The temperature related parameter has the fixed value s = 0.1. 
Now, the reset state is not the pure state Preset = | + +}{++ but the mixed 
state preset = p\ + +) (+ + | + (1 — p)/21. The three curves, from top to bottom, 
correspond to p = 1, p = 0.98, and p = 0.97. With a mixed reset state, a third 
threshold for large but finite reset rate r appears above which the state is not 
entangled in contrast to the case of pure reset states where the entanglement 
vanished only for r — > oo (topmost curve). 



The picture that emerges from all these results is the following. Entanglement can 
prevail in dissipative, open quantum systems that are far away from thermodynamic 
equilibrium. For gas-type systems treated in this section, a reset mechanism can evoke 
steady-state entanglement even for infinite temperature of the environment generically, 
i.e., independently of the specific form of the interaction Hamiltonian or decoherence 
channel. 

In the next section we show that steady-state entanglement appears also in 
strongly coupled systems with an appropriate reset mechanism. 



3. Strongly coupled systems 

In gas-type systems, we can treat the local noise channels separately for each qubit 
as explained above. If these local channels correspond to a heat bath, they drive 
each qubit individually to the thermal state of the local, free Hamiltonian, i.e., they 
populate the eigenstates of the free Hamiltonian according to the Boltzmann factor. 
Although the effective interaction Hamiltonian in the master equations is represented 
as continuously acting, the physical interaction process in gas-type scenarios is viewed 
as a short collision event. Hence, the interaction Hamiltonian does not influence the 
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energy spectrum considerably and does not modify this thermal state. In strongly 
coupled systems, on the other hand, interactions of quanta of the heat bath with the 
system qubits affect the system as a whole. In this sense, the decoherence process acts 
globally on the system, inducing transitions between joint eigenstates. In this section 
we will shortly discuss the master equation describing a strongly coupled spin system 
in contact with a thermal, photonic bath. We will see that the resulting equilibrium 
state, the thermal state, can be entangled below a certain temperature threshold if 
the ground state of the Hamiltonian is entangled. When we add a reset mechanism we 
find that, in contrast to the gas-type scenario, entanglement in the steady state can 
exist only below a certain temperature threshold. However, the novel feature is that 
this threshold is typically much higher than for the thermal state. Finally, we describe 
the influence of the master equation parameters on the respective steady state, and, 
with this insight, formulate a general condition under which a reset mechanism can 
lead to an entangled steady state. 

Let | a) (| b)) be momentary eigenstates of some non-degenerate system 
Hamiltonian H(t) with eigenenergies ui a (wfejjl We define N a b :— ^e^( Ua ~ Ui ^ — 1) _1 
with [3 = 1/T being the inverse temperature. Often N a b is written as n, and we 
explained the connection to the parameter s in the last section. The master equation 
for a spin system, coupled with strength 7 to a heat bath consisting of photons, is [12] 

p = - T ^[Ar &a |g ba | 2 |< a |^)|6)| 2 + (A^ afo + ^l^pKfol^la)! 2 ] 

j, a,, b 

x{\a)(a\p + p\a)(a\-2(a\p\a)\b)(b\}, (12) 

Here, g{oj) is the spectral density, for which gt, a — g{^b — u a ) if uib > oj a and gb a = 
else. For small system one can justify to treat the spectral density as constant {g = 1 
if iOb > u a ) and merely tune the overall coupling constant 7 [13) . Observe that we did 
not include non-radiative contributions as opposed to the master equation ([3]) with 
noise terms (j4]). 

The master equation (|12|) drives any initial density matrix to the thermal state 
of inverse temperature (3. That means, the ground state and also the excited states 
are populated according to the canonical distribution. 

We will study the master equation (|12[) for an Ising Hamiltonian with transverse 
magnetic field, briefly discuss the solution without reset mechanism (thermal state), 
and then turn to an analysis of the full master equation with reset mechanism. 

3.1. Master equation without reset 

An Ising Hamiltonian with transverse magnetic field has the form 

H^giaP^+bia^+a^)}, (13) 

and the eigenvalues are — gVl + 4b 2 ,—g,g,g\/l + Ab 2 with corresponding eigenvectors, 
expressed in the standard basis, 

|Vo) = N(l, (-1 + v / l + 46 2 )/26, (-1 - v / l + 46 2 )/26, 1), 

|Vi> = 1/72(0,-1,1,0), 

IV2) = l/\/2(-l, 0,0,1), 



X Note that h = ks = 1 throughout the paper. 
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|^3> = JV(1, (-1 - v / l + 46 2 )/26, (-1 - v / l + 46 2 )/26, 1) 



(14) 



where N = (2 + l/2|(— l + yT+ 46 2 )/&| 2 ) -1 / 2 provides normalization. We exclude the 
case 6 = where the ground state would be degenerate, a case not properly described 
by the master equation (fT2|) . The ground state of this system is the first eigenvector 
in equation (fT4"|) . Since this state is entangled, so is the thermal state below a certain 
temperature threshold. We see this directly from the negativity of the thermal state 
Pthcrmai = cxp{-/3.ff }/ tr exp{-/3.ff} , which is 



For any fixed g and b, the non-trivial part in this formula goes to the value — 1/4 when 
(3 — > 0, while it goes to l/(2%/46 2 + 1) > for /3 — > oo. The threshold value for /?, 
where the non-trivial part becomes exactly zero, can be easily computed numerically 
for any given parameters g and b. In terms of /3, the thermal mixture of the eigenstates 
is separable below this threshold, which, in terms of T — /3 _1 , means that the mixture 
is separable above that critical temperature. From [15] one can see that the critical 
temperature grows linearly with g and monotonously, but sub-linearly with b. 

3.2. Master equation with reset 

We keep the Ising Hamiltonian with transverse magnetic field as above, but extend 
the master equation ([12"]) by the reset term 

N 



with N — 2 and reset state = |+) for both qubits. We solve the resulting master 
equation numerically. In Figure[5]we see how the entanglement for different reset rates 
r develops over time t from the value zero in the initial product state p = | + +)(+ + | 
to its final value in the steady state while all other parameters are kept fixed (see 
figure caption) . We notice that small, non-zero reset rates decrease the entanglement 
in the steady state until it is gone, while larger rates can bring entanglement back. 

3.2.1. Influence of the parameter r. To explain this effect, recall that the master 
equation mimics the averaged density matrices that would be obtained from (infinitely 
many) simulation runs of the system. The reset rates of the master equation are 
related to probabilities that in a simulation a reset took place during a certain time 
interval. Although our reset processes are strictly speaking local, let us assume for 
the sake of argumentation that the reset happens on both qubits simultaneously, thus 
effectively restarting the process again from the beginning whenever a reset occurs in 
a simulation. For small rates r, i.e., for small probabilities that a reset takes place in 
the simulation, the system can come close to its thermal equilibrium state before it is 
reset. When we average the density matrices over many simulation runs, we average 
matrices that are mostly close to the unique thermal equilibrium state, and hence 
also the mixture will still retain entanglement. When the rates get larger, the density 
matrices over which we average become more and more diverse since they will be far 
from equilibrium and fluctuations occur. As a consequence the average density matrix 
will have no entanglement. When the reset rate is above a certain threshold, we will 




(15) 




(16) 
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find entanglement in the system again (as we did in the gas-type systems) because 
now the density matrices over which we average become similar again. Now, they are 
close to the state that has unitarily evolved for a time of order 1/r from the initial 
reset state. In the limit r — > oo the state is constantly kept in the initial product state 
with zero entanglement. In this way we can understand how the two entangled regions 
arise. The first is an artifact of the entangled thermal state that is more and more 
destroyed by the reset mechanism. This first region could also be present in a gas-type 
model. The second region is the one that is really created by the reset mechanism 
just as in the gas- type model. We can directly see in Figure [3] that entanglement in a 
thermal state is a truly different effect from entanglement that is created by the reset 
mechanism. 



(a) (6) 




Figure 9. (a) Solution of the master equation l|12[ l with reset term (|16fl for a 
strongly coupled system. The unit timescale is 7 -1 and the time t is plotted in 
these units. The temperature parameter was chosen as /3 = 1000, and, at this 
low temperature, the steady state of J 1 2H is entangled even without reset, the 
other parameters being g = IO7 and b = 0.1. When increasing the parameter r 
of the reset mechanism, the steady-state density matrix changes, as explained in 
the text, fist becoming separable and then entangled again, (b) Cut through the 
same plot for different r. 



3.2.2. Influence of the parameter 7. Although the two effects are truly different, this 
does not mean that for certain parameter regimes, the two regions cannot overlap, 
see Figure ITOl Imagine the coupling to the photon bath, 7, is increased. This means 
that the system is driven towards thermal equilibrium faster than before. Hence, 
following the arguments from above, the system can tolerate higher reset rates before 
the entanglement in the first region is destroyed. The stronger coupling to the photon 
bath suppresses the entanglement in the second region, and as an overall effect we 
see that the two regions need not be separate. Note that the entanglement in the 
thermal state for r = is independent of 7 because then it does not matter how fast 
equilibrium was approached. 

3.2.3. Influence of the parameter b. The transverse magnetic field with relative 
strength b splits up the energy levels of the Hamiltonian (we exclude the degenerate 
point 6 = 0). The ground state will contain less and less entanglement as b increases 
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Figure 10. (a) Influence of the parameter 7, which describes the strength of 
the coupling to the photon bath. Here, we choose g~ 1 as unit time. The other 
parameters are j3 = 10 and b = 0.1. (b) Cut through the same plot for different 
values of 7. 

(and will approach the product state | ) for b — ► 00). Hence, at zero temperature, 

the entanglement in the equilibrium state for r = will go down for increasing b (see 
red line in Figure UlTa)). For a thermal state with T ^ 0, i.e., (3 ^ 00, the ground 
and exited states get mixed. When b is small, the splitting between ground and first 
exited state is small, and the mixture will be close to a separable state. For increasing 
b, the larger energy split leads to an increased population of the ground state relative 
to the exited states at the same temperature 1//3, and the entanglement will increase. 
When b gets even larger, the thermal state will be close to the ground state, but we 
know, that the ground state for large b is only weakly entangled. Hence, there will 
be some b for which the entanglement in the thermal state is maximal (see red line 
in Figure [TTTb), and the identical line in Figure [TTTc)). When we switch on the reset 
mechanism, we see that an increasing r destroys the entanglement in the thermal state 
as before, but for increasing b the regions one and two (artifact of thermal state vs. 
true reset state entanglement) quickly overlap. That is because the magnetic field 

tends to drive the state towards | ), whereas an increasing reset mechanism drives 

the state towards | + +)• If we choose the reset state as |— ), the two effects do not 
compete and both drive the state towards an unentangled product state (see light grey 
areas in Figure fTTTc) and compare to the light grey areas in Figure [TTTb)). How fast 
an increasing r destroys the entanglement in the thermal state is almost independent 
of the entanglement in the thermal state. For larger r the influence of b plays less and 
less a role and there is almost no difference between Figures [TlTb).(c') for large r. 

3.2.4- Influence of the parameters g and (3. The coupling strength of the Hamiltonian 
g also splits the energy levels. Hence with increasing g there will be more entanglement 
in the thermal state at some finite temperature. Again, the speed with which 
an increasing reset rate destroys the entanglement in the thermal state is almost 
independent of g (see Figure HUJa)). Most interesting is the influence of the 
temperature 1/(3. Figures [T27a)-(c) show plots which contain information about the 
equilibrium-state entanglement for different temperatures (r = 0). We see how the 
thermal states get less and less entangled for increasing temperatures, so that region 
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Figure 11. Influence of the relative magnetic field b. The unit timescale is given 
by 7 _1 and g = 5O7. Plots (a) and (b) show the difference between (almost) zero 
temperature (/3 = 10000) and some finite temperature ({3 = 0.2). The behavior of 
the red curves, for r = 0, are explained in the text. Plots (b) and (c) demonstrate 
the influence of different reset states in connection with the magnetic field b. Plot 
(b) has reset state |+), (c) has reset state |— ). 



one vanishes quickly as expected. But, for the same temperatures, the reset rate r 
can still produce entanglement in a steady state! On the other hand, for fixed g, there 
is some temperature threshold above which no reset rate can produce entanglement. 
This threshold is in contrast to the case of gas-type systems. 

(o) (6) (c) 

f3 = 1000 (3 = 0.2 (3 = 0.02 




Figure 12. Negativities for increasing temperature (decreasing /3 = 1/T). The 
unit timescale is given by 7 , the relative magnetic field by b = 0.1. The 
parameter region in the r — g plane, where steady-state entanglement occurs, 
becomes smaller for increasing temperatures. However, the temperatures, for 
which entanglement can exist with reset mechanism, are much higher than the 
temperatures, for which the thermal state is entangled. The red line in (c) 
corresponds to the red line in Figure 1131 (see that figure caption and the text) 
and is drawn for comparison only. 
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3.3. General conditions for steady-state entanglement 

This discrepancy between gas-type and strongly interacting systems raises the deeper 
question: What are the conditions under which the reset mechanism can create 
entanglement in the steady state? We already see that the reset mechanism is different 
from a cold bath - equivalently the replacing of system particles by fresh, standard 
ones - because we would expect that some cold bath can always counteract the 
influence of the hot photon bath. The condition that the solution of the master 
equation (a completely positive (CP) map) at r = is entangling at some point in 
time is certainly necessary, since, as pointed out before, the reset mechanism does not 
introduce entanglement itself. The question is: Is a solution of the master equation 
at r = that creates entanglement on some short time scale sufficient such that the 
solution for some r > is a CP-map with entangled steady state? Unfortunately this 
is not true. The reset rate r itself can influence the solution of the master equation 
in such a way that although the solution for r = was entangling on some short time 
scale, the solution for larger r need not be. 

The condition for steady-state entanglement is the following. If the solution of 
the master equation for some r > is entangling on a time scale of order 1/r then 
the steady state of this solution will also be entangled. For illustration, we think once 
more in terms of many hypothetical simulation runs. As stated above, the states over 
which we have to average will be close to the state that has unitarily evolved for a time 
1/r from the initial reset state. When this state has a certain amount of entanglement, 
then so does the mixture of states close to it (see Figure [13]) . 



Figure 13. Illustration of the condition under which a reset mechanism can create 
entanglement in the steady state. At time -ft = 2 the state is already close to the 
true steady state, so the red line corresponds to the red line in Figure [I2{c) since 
g = IO7 in this plot, and the other parameters are the same as in Figure IT^ c) . 
The reset can create entanglement in the steady state if, for given r, there is 
entanglement in the time-evolved state at time t cx 1/r (see text for details) . The 
green curve is the curve t = 2/r illustrating this result. 
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The condition does not only hold for the specific Hamiltonian or the specific heat 
bath chosen here. It is valid for a large class of Hamiltonians (with appropriately 
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chosen reset states) and baths. 

In the next section we treat the multipartite case. 

4. Multipartite case 

For multipartite spin systems, we will show that the reset mechanism can create 
steady-state entanglement in a similar way. The parameter regions where this happens 
are comparable to the 2-qubit case. The values of the negativity, or rather its 
generalization, the average negativity, stay almost constant with increasing system 
size. Note, however, that larger systems could have larger negativities, so if we divide 
the actual (constant) negativity by the maximal possible negativity, then this quantity 
would go down for growing system size. We will also consider entanglement in reduced 
density matrices. Since in a reduction from ./V to, say, 2 qubits the traced out N — 2 
qubits act as an additional noise source, it is not surprising that the parameter region 
where we find steady-state entanglement in the reduced systems shrinks with growing 
N. If the number of particles fluctuates according to some distribution, our best 
description of a reduced density matrix is a mixture of reductions originating from 
different system sizes. It is remarkable that even in this case there is some parameter 
region, where steady-state entanglement is found. 

In the following, we motivate and explain the entanglement measures we are 
using and then demonstrate the above features in both gas-type systems and strongly 
interacting systems. 

4--1- Gas-type systems 

It is straightforward to generalize the master equation for the gas-type system, 
equation ([7} to N > 2 qubits. The relevance of the entanglement quantities we 
are going to use needs to be motivated, though. We turn once more to the example 
of a spin gas. Imagine that the spin gas is in a box of volume V. This box has one 
semi-permeable wall through which particles can leave the box, and another through 
which particles from a "cold reservoir" can enter. By "cold reservoir" we mean that 
the quantum state of the particles in this reservoir is in some sufficiently pure standard 
state. The motional degrees of freedom of theses particles, on the other hand, are in 
thermal equilibrium with the outside environment just like the system particles in 
the box. Assume that the density of the gas in the box is n. Then there are on 
average r/V =: A particles in the box. The distribution of the number of particles 
that are in the box is a Poissonian p\(n) = e~ x \ n /nl. When we observe the spin gas 
after certain time intervals, which should be long enough such that the gas always 
reaches its equilibrium state, we sample the distribution and get information about 
the density matrices with a corresponding number n of qubits. The density matrices 
we can reconstruct after we collected a certain, sufficient amount of information is 
close to the steady-state density matrix of the master equation for n qubits. We 
are interested in the entanglement properties of the gas and we will look at different 
aspects of entanglement in the following. 

All these aspects of entanglement are quantified by measures that are based on 
the negativity or the average negativity. The average negativity Af is the negativity 
averaged over all possible bipartitions of the system [9] . Non-zero average negativity 
ensures the presence of some form of entanglement in the system. 

Specifically, we study three types of entanglement: 
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(i) The average negativity of the density matrices with n qubits, averaged over the 
Poissonian distribution of the number of particles in the system: (J\(p n ))px(n)- 

(ii) The negativity of reduced 2-qubit density matrices averaged over a renormalized, 
truncated Poissonian distribution: (J\f(p n _ »2))p A (n>2)- 

(iii) The negativity of averaged, reduced 2-qubit density matrices: Af(p). 

Now, we lay out, what these quantities mean, and which aspect of entanglement 
they describe. 

(i) If we ask how much entanglement we find in the system on average we are led 
to the quantity (Jv(pn)) PX (n)-> which is the expectation value of the average negativity 
of density matrices with different n, where p\(n) is the Poissonian probability 
distribution for the number n of particles in the system. Observe that if we 
disallowed the fluctuation of the particle number in the system and introduced the 
reset mechanism by other means (e.g. measurement and decay to standard state inside 
the box), the quantity of interest would simply be M{pn) for fixed system size N. 

(ii) When we look at subsystems we are led to slightly different quantities. Let us 
fix the subsystem size to 2 qubits. We call the reduced density matrices of originally n 
qubits pn~,i where we assume n > 2. In gas-type systems there can also be zero or one 
particle in the box (especially if 77 cx A is small), and the entanglement is simply zero 
in these cases. Since a "reduction" of a one or zero-qubit density matrix to a 2-qubit 
density matrix makes no sense, we simply exclude these cases and rescale the truncated 
Poissonian p\ [n > 2) to the distribution p\(n > 2). The quantity (J\f(p n ->2))p x (n>2) 
therefore tells us how much entanglement a subsystem of two qubits contains on 
average (for 2 qubits N = AT) . 

(iii) When we look only at a subsystem of two qubits disregarding the number of 
particles n in the system, then our best description of the 2-qubit density matrix is 
the average density matrix p := (p n -*2)p x (n>2) with entanglement N{p). 

When we compare the three kinds of entanglement defined above we see that 
the conditions that one of these quantities be non-zero are increasingly stringent. To 
find entanglement in the reduced system is a more stringent condition since tracing 
out the other particles has the same effect as an additional noise source. Also, to 
find entanglement in the averaged density matrix p is a stricter condition since the 
averaging increases entropy, i.e., tends to make the matrix more mixed. If we keep 
the particle number fixed, N > 2, there is just the quantity Af(pN->2) to describe the 
entanglement in the reduced state. 

To compute the average negativity is a hard task. The system size, and hence the 
number of differential equations we must solve, and also the number of bipartitions 
scale exponentially. To simplify the computation, we consider a symmetric situation, 
where all qubits interact pairwise via Ising interactions and are subject to dephasing 
noise (J5)). In Figure[l4]we plot the three entanglement measures (N{pn))p x {n) (blue, 
dashed curve), (J\f(p n -+2))p x ( n >2) (orange, solid curve), and A/"(p) (black, dashed- 
dotted curve) for A = 2 and for a qubit number that fluctuates between two and five. 
Since the meaning of these measures is different, one cannot compare the absolute 
values represented by the curves directly with one exception. The points, where the 
curves become non-zero, must, from left to right, appear in the order explained in the 
previous paragraph, i.e., blue first, representing measure (i), orange next, representing 
(ii), and black last, representing (iii). 
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Figure 14. Measures (Af(pn,)) Px („) (blue, dashed), (M{p n -,2))p x (n>2) (orange, 
solid), and Af(p) (black, dashed-dotted) as functions of the reset rate r/7. 
The particle number fluctuates between two and five, and the fluctuations are 
accounted for by a Poissonian or truncated Poissonian weighting with A = 2 as 
explained in the text. The interaction strength of the Ising interactions between 
all particles is g = 2O7, while the strength of the local, free Hamiltonian |[T) is 
u> = 5O7. Kinks in the orange curve stem from the averaging over negativities 
whith different supports (see (ii) in the text). 



4-2. Strongly coupled systems 

Eventually, we study equations (fT2"|) and (fTB|) in the multipartite case. We obtain 
similar results as in the case of gas-type systems underlining again how generic the 
reset mechanism is. 

Although a fluctuation of particles may seem less natural in strongly coupled 
systems as compared to gas- type systems, we will look at the exact same entanglement 
quantities, so that the corresponding plots are directly related to each other. As 
Hamiltonian, we choose a sum of pairwise Ising interactions and magnetic fields in x 
and z direction according to 



where the small gradient magnetic field in z direction is introduced for technical 
reasons to lift degeneracies in the Hamiltonian. Figure [15] shows a plot of the same 
entanglement measures as in the previous subsection. As in the gas-type scenario, 
the feature that entanglement can be created by a reset mechanism holds also in the 
multipartite case. While we have shown the genericity of the reset mechanism with 
respect to the Hamiltonian and noise process already in previous sections, here we 
demonstrate that the reset mechanism is also generic with respect to system size. 

5. Summary 

We have shown that entanglement can be present in dissipative, open quantum systems 
far from thermodynamic equilibrium if we assume the existence of an additional 




(17) 
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Figure 15. Measures (Af(pn)) Px ( n ) (blue, dashed), (A/"(pn->2))p A (n>2) (orange, 
solid), and Af(p) (black, dashed-dotted) as functions of the reset rate r/7. 
The particle number fluctuates between two and five, and the fluctuations are 
accounted for by a Poissonian or truncated Poissonian weighting with A = 2. The 
other parameters of the Hamiltonian II17I I are g = 157 and b = 0.1. The inverse 
temperature is /3 = 0.2. Kinks in the orange curve stem from the averaging over 
negativities whith different supports (see (ii) in the text). 



mechanism that "resets" the particles, at a certain rate, into a single-particle, low- 
entropy state. For a 2-qubit toy model of a gas-type system, we have analytically 
solved the master equation consisting of a Hamiltonian part, a noise channel, and 
the proposed reset mechanism. For special cases we have been able to give closed 
expressions for the entanglement as a function of the parameters of the master 
equation. We have extended the analysis to similar models with other interaction 
Hamiltonians, decoherence models, and imperfect reset mechanisms. We have treated 
the situation of strongly correlated systems by the same means and we have given 
conditions under which steady-state entanglement arises in this case. Finally, we have 
shown that in systems consisting of more qubits, and even in systems with fluctuating 
particle number, steady-state entanglement can prevail. 

Many systems are conceivable for an experimental realization of such dissipative, 
open quantum systems. For instance, one may consider ions in microtraps that interact 
via an induced dipole moment |14j leading effectively to a continuously operating Ising 
interaction. Decoherence, dominated by dephasing noise, appears naturally in such 
systems, and the reset mechanism may, e.g., be achieved by periodically applying a n- 
pulse that couples the internal level |1) to a metastable auxiliary level \a) that decays 
rapidly to |0). The state afterwards is always |0), which can be mapped to |+) by a 
subsequent Hadamard operation. 

For charge manipulated quantum dots, the exchange interaction leads to a 
continuously operating Heisenberg interaction between neighboring electron spins by 
lowering the potential barrier [15] . The effect of surrounding nuclear spins may be 
described by dephasing noise, while the reset mechanism can consist in replacing an 
electron by a fresh one from the surrounding Fermi sea, prepared in a suitable state 



(e.g., |0». 
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Atomic beams interacting via a cavity mode |16| may also serve as a toy example 
of such systems far from thermodynamic equilibrium. 

Note that a reset mechanisms could be realized in many physical ways, including 
a measurement with subsequent preparation of the state, coupling or decay to 
metastable auxiliary states, as well as replacing a qubit by a fresh one. While 
these suggestions for implementations aim at demonstrating how a reset mechanism 
would be realized in experiments, the effect itself is generic and other realizations 
are conceivable. In particular, one might try to find such a reset mechanism in less 
controlled, maybe even biomolecular systems consisting of many particles. 
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Appendix A. 

In this appendix we derive the solution of the master equation for the Ising 
Hamiltonian ^ explicitly. In the final formulae we restrict ourselves again to the 
dephasing channel ([5]). 

Instead of solving the master equation by spectral decomposition of the matrix 
A associated with the Liouville operator £ (see section [2~4l) . we solve the system of 
linear differential equations step by step. Since not all of the differential equations are 
coupled, this leads to simpler expressions. The disadvantage is, however, that the set 
of solutions contains 9 integration constants that must be determined afterwards for 
given initial conditions. 

We expand the 2-qubit density matrix p in the standard basis |s) = |siS2) with 
Sj e 0,1 and <r£ |s) = {— l) Sj |s) . The expression for the density matrix becomes 
p = J2s' 1 s' 2 s 1 s 2 Cs' 1 s' 2 sis2 l s i s 2)( s i s 2|- Inserting this expansion into the master equation, 
and defining the two functions 



Here, the operation © means addition modulo 2. Fortunately, these 16 differential 
equations are not fully coupled. The coefficients Cooocb Coioi> Cioio> an d Cim on the 
diagonal of the density matrix are coupled only to themselves. Once we have solved 
these equations, we can treat the diagonal coefficients as known inhomogeneities in 
the other equations. The off-diagonal coefficients Cqooi and Cion are coupled among 




- (x®l)), 



d t C s , s , SlS2 = {-ig (-1)^+4 _ 




+ AOi) si) + M4, s 2 ) ~ 2r}C s;4siS2 

+ /2(Sl> s l)^(siel)s 2 (siel)s 2 + /2( s 2> s 2)C s ' i ( s ^i) Sl ( S2e i) 
+ r /2{Cos 2 Os 2 + Cls 2 ls 2 + CsiOsiO + C^i^i}. (A.l) 
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themselves and to the diagonal, so we can solve them next. The same is true for the 
pair Cooio and Com- Finally, the anti-diagonal coefficients Coon and Cono are cou- 
pled to Coooi,C'ioii,Cooio: and Com, or to their complex conjugates. We solve these 
as a last step, and all other coefficients are given by the Hermiticity of the density 
matrix. The solution is now straightforward in principle. However, the expressions 
for the matrix coefficients are still space-consuming, so we will give them only for the 
special case of the dephasing channel. 

For the dephasing channel ([5]), the structure of the differential equations is still 
the same and not simplified, we save space only because we have fewer parameters 
and a symmetric situation. The solution for the diagonal elements then reads: 

Coooo = \ + \D2e- 2rt + \Dze- rt 
Coioi = \~ \D 2 e- 2rt + \Oie- rt 
C 10 io = \~ \D 2 e- 2rt - \v±e~ rt 
C lm = \ + \D 2 e- 2rt - \D^ rt 

The integration constants D 2 , D3, and D4 accommodate the initial conditions. 
Since trp = 1 is a constraint, there is no free constant D\. 
The off-diagonal elements are: 

r(-ig + r + j+f ) (D 3 + D 4 ) e - rt r(2ig - 2 7 - iu)) 

L'OOOl 



— e 

X 



4 (2g 2 + (r + 7 + f ) (r + 2 7 + iu))) 4 (4<? 2 + (2 7 + iuj)(r + 2 7 + iu;)) 

jt(V+47+2iw+<y/r 2 -16g 2 ) 

f e V r2 - 16 9 H OB 2 i^r 2 - 16g 2 - Aig^j + ODi Uig + ^r 2 - 16 5 2 ) 



„ _ r(ig + r + 1 +f) (D 3 + D 4 )e- rt {2ig + 2j + iuj)r 
< loii — " _ ~ — 



4 (2g 2 + (r + 7 + f) (r + 2 7 + iu))) 4 (4g 2 + (2 7 + iu))(r + 2 7 + iu))) 

e -|t(3r+4 7 +i( 2w +Vl6^-H)) r / QDi _ el VT%^t D 2 



The coefficients Cooio, Com are very similar to C0001, Cion, except that D4 must 
be replaced by -D4, and the integration constants are OD3, OD4 instead of ODi, OD2. 
Finally, the elements on the anti-diagonal are 

n = e -2t(r+27+M An ■ (2r + 2j + iu))r 2 

0011 1 4 (4g 2 + 2r 2 + (2 7 + iu;) 2 + r(6 7 + 3m;)) (r + 2 7 + iu) 

ie~ rt gD 3 r 2 ^^t(^3r+4i+2iiu+iy/l6g 2 ~r 2 ) 

(4g 2 + (27 + iu>)(r + 2j + iu)) (r + 4j + 2iu) 



x r 



4.g + ir + 0V - r 2 (ODi + OD 3 ) 



+ 4i 7 - 2u> + \/l6g 2 - r 2 
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e W^9 2 -r 2 t (_ Ag _ ir+ y/Tfyf - r A (0D 2 + 0D 4 )\ 



— ir — 4ij + 2oj + y^lQg 2 — r 2 I 

and 

Coon = e t( -- 2r -^AD 2 

r 2 (8(r + 7)5 2 + (r + 2 7 ) (4(r + 7) 2 + oj 2 )) 
4(r + 2 7 ) (w 4 + (-8.g 2 + 5r 2 + 8 7 2 + 12r 7 ) us 2 + 4 (2g 2 + (r + 7 )(r + 2 7 )) 2 ) 
ie~ rt g (i(r + 4j)ujD 3 + (4g 2 - uj 2 + 2 7 (r + 2 7 )) D 4 ) r 2 



(r + 4 7 ) (16.g 4 + 8 (2 7 (r + 2 7 ) - w 2 ) .g 2 + (4 7 2 + lu 2 ) ((r + 2 7 ) 2 + w 2 )) 
i / / — / f 4g - ir + \/l6g 2 - r 2 ) ODi 

e -it(^3r+47-2iw-i\/16g 2 -r 2 j^/ V V / 

\ — — 4i7 + 2lj + yflQg 2 — r 2 

e -*V le 9 2 -r 2 t ( Ag + ir+ ^ lgg 2 _ r 2^ QD 2 \ 

ir + 4i7 — 2lu + y^ldg 2 — r 2 I 
, c / — ^-n / f 4g + + a/16o 2 - r 2 ) OD 3 

e -it(3r+47+2iw+i % /l6g 2 -r 2 J^/ V y ' V y 

I ir + 4«7 + 2w + \A6ff 2 - r 2 



eV 16 f 2 - r2 * f-4g - ir + ^16.g 2 - r 2 ) OD 4 \ 
— ir — 4i7 — 2lo + \/ 16g 2 — r 2 / 

All other coefficients follow from the Hermiticity of the density matrix. The form of all 
matrix coefficients is similar. First, there are the parts with the integration constants 
that fall off exponentially with time. The characteristic exponents are the eigenvalues 
of the homogeneous parts of each linear differential equation or system of equations 
(multiplied with time). As pointed out in 12.4. H these exponents are the spectrum of 
the total Liouville super-operator defined by p = Cp with values 

{0, -r, -2r, -2(r + 2 7 ), -2(r + 2 7 ± iu), 

-l/2(3r + 4 7 + y/-16g 2 + r 2 ± 2iu>, -l/2(3r + 4 7 - ^/ -IQg 2 + r 2 ±2iuo} 

and multiplicities {1, 2,1,2,1 + 1,2 + 2,2 + 2} respectively. Second, there always is a 
part independent of t (belonging to the eigenvalue 0) that represents the value in the 
steady state. 

Appendix B. 

In this appendix we present the steady-state solution for the master equation ([7J with 
Ising Hamiltonian, local noise channels as in equation (01 and reset states = |+) 
to |Appendix B[ The solution in form of the matrix coefficients in the computational 
basis is given by: 

_(r + 2Bs) 2 

L'OOOO 



C0101 — C1010 — 



4(£ + r) 2 

(r + 25(1 - s))(r + 2Bs) 



A(B + r) 2 
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_ (r + 2i?(l-.s)) 2 
Clm " 4(B + r)* 

_ r(r + 2Bs){B + C - 2ig + 2r + iw) 

0001 ~~ A{B + r) (V + (C + r + iuj){C + 2r + iw) + B{C + r + 2ig(2s - 1) + ico)) 

Coon = {4(B + r) (4g 2 + (C + r + iu)(C + 2r + iw) + £(C + r + 2£ 9 (2s - 1) + ioj)) 
x(C + r + iuj)} |r 2 (B 2 + (C + 2i 5 + 3r - 4i#s + iw)S + r(C + 2r + iuj) \ 
Cono = {4(B + r){C + r) (uj 4 + (2C* 2 + 6rC - &g 2 + 5r 2 ) uj 2 
+ (4g 2 + (C + r)(C + 2r)f + 2B ((C + 2r)uj 2 + 2g(2C + 3r)(2s - l)w 
+ {C + r) (4g 2 + (C + r)(C + 2r))) + B 2 ((C + r) 2 + ( 5 (4s - 2) + w) 2 )) 
x |r 2 ((C + r)S 3 + ((C + r)(2C + 5r) - 16g 2 (s - l)s) B 2 

+ (C* 3 + 7rC 2 + 4.g 2 C + 14r 2 C + 8r 3 + (C* + r)w 2 + 12.g 2 r - Ag{C + r)(2s - l)w) B 
+r(C + r)w 2 + r(C + 2r) (4.g 2 + (C + r)(C + 2r))) }. 

Furthermore, Cooio = Coooi: and Com — Ciou are obtained from Coooi by replacing 
g — > — gr, s — » (1 — s) in the numerator. The other coefficients are given by the 
Hermiticity of the density matrix. Observe that the dephasing and the depolarizing 
channels are included as special instances of the parameters B,C, s in this analytic 
expression. The dephasing channel results from putting B = 0, s drops out, and 
renaming C = 2-f. The depolarizing channel is given by s = 1/2, B = C, and 
renaming C = 47/3. 

The plot in Figure 0] is based on this solution. 
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